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Abstract  Under  conditions  of  isostaticity  in  granular  media, 
the  contact  forces  for  all  particles  are  statically  determinate 
and  forces  can  be  computed  without  recourse  to  deformation 
equations  or  constitutive  relationships.  Given  that  stresses 
represent  spatial  averages  of  inter-particle  forces,  the  stress- 
equilibrium  equations  for  the  isostatic  state  form  a  hyper¬ 
bolic  system  of  partial  differential  equations  that  describe 
the  internal  stress  state  using  only  boundary  tractions.  In  this 
paper,  we  consider  a  Cosserat  medium  and  propose  closure 
relationships  in  terms  of  stresses  and  couple  stresses  from 
observations  of  stress  variations  in  the  critical  state  regime 
from  discrete  element  simulations  and  experiments  on  sand, 
even  though  the  isostatic  condition  is  only  satisfied  in  an 
average  sense.  It  is  shown  that  the  governing  equations  are 
hyperbolic,  which  can  be  solved  using  the  method  of  char¬ 
acteristics.  Examples  of  both  analytic  and  numerical  solu¬ 
tions  are  presented.  These  examples  clearly  demonstrate  that 
stress  chains  (characteristic  lines)  form  oblique  angles  with 
the  assumed  direction  of  the  force  chains. 
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1  Introduction 

It  is  well  known  that  in  two  dimensions  the  equations  for 
classical  plasticity  form  a  hyperbolic  system  of  partial  dif- 
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ferential  equations  (e.g.  [1]).  These  equations  are  obtained  by 
rendering  the  system  of  momentum  balance  equations  stati¬ 
cally  determinate  via  the  yield  criterion  (e.g.  Coulomb-Mohr 
relation  in  [1]),  which  provides  an  additional  relationship 
among  stresses  necessary  to  close  the  system  of  equations. 
Thus,  the  stress  in  a  medium  undergoing  plastic  flow  can  be 
computed  from  the  boundary  conditions,  without  recourse  to 
the  deformation  or  stress-strain  relationships. 

Blumenfeld  and  co-workers  [2-5]  examined  isostaticity, 
a  state  in  which  the  forces  are  statically  determinate  for  an 
assembly  of  rigid  particles.  Continuum  equations  derived 
from  this  state  reflect  the  conditions  at  the  particle  scale  and 
are  likewise  statically  determinate,  giving  rise  to  a  set  of 
governing  equations  that  are  hyperbolic.  In  reality,  all  parti¬ 
cles  possess  some  compliance,  form  frictional  contacts  and, 
for  dense  granular  materials,  are  generally  hyperstatic  in  the 
stable  pre-failure  regime.  The  question  that  then  arises  is 
under  what  states  might  the  condition  of  statically  determi¬ 
nate  stresses  or  isostaticity  be  reasonably  assumed  in  a  con¬ 
tinuum  as  well  as  at  the  particle  scale?  Before  proceeding, 
note  that  the  issue  of  compliance  for  frictionless  (Hertzian) 
contact  was  addressed  by  Blumenfeld  [6] ;  therein  he  showed 
that  provided  the  mean  coordination  number  in  d-dimensions 
is  d  +  1,  a  macroscopic  system  of  compliant  but  frictionless 
particles  can  be  mapped  onto  an  equivalent  assembly  of  per¬ 
fectly  rigid  particles  that  approximately  supports  the  same 
stress  field.  The  error  or  discrepancy  between  the  two  stress 
fields  decays  with  system  size  according  to  a  power  law. 
In  this  paper,  however,  we  focus  on  developing  a  statically 
determinate  set  of  Cosserat  constitutive  relations  for  a  gran¬ 
ular  assembly  comprising  frictional,  compliant  particles. 

It  is  instructive  to  explore  the  differences  and  similarities 
between  the  system  of  equations  from  classical  plasticity  and 
those  of  isostaticity  in  [2-5].  As  in  two-dimensional  plastic¬ 
ity  theory  (e.g.  [1]),  the  stress  in  an  isostatic  medium  can 
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Fig.  1  Stress  chains,  force  chains  (dark  gray  particles ),  and  the  rela¬ 
tionship  between  stress  along  a  characteristic  line  and  the  average  stress 
of  the  sample  volume  (left  and  top  right).  Depiction  of  a  growth  and 
collapse  by  buckling  of  force  chains  in  a  shear  band  (bottom  right): 
a  force  chain  grows  in  upper  band  region  from  a  particle  chain  under 
compression  opposing  the  shear,  while  an  existing  force  chain  in  lower 
band  region  undergoes  buckling  with  middle  segment  under  extension 
and  rotating  in  the  direction  of  shear 


be  found  solely  from  the  tractions  applied  to  the  boundaries. 
Thus  the  key  difference  between  the  equations  of  plasticity 
and  isostaticity  lies  in  the  fact  that  while  the  former  applies 
only  when  the  continuum  is  in  an  ideally  plastic  state,  the  lat¬ 
ter  applies  at  all  stress  states  and  for  both  continuum  and  par¬ 
ticle  scales.  Another  aspect  that  is  common  to  both  systems  of 
equations  is  their  type,  i.e.  hyperbolic.  The  physical  interpre¬ 
tation  of  hyperbolic  partial  differential  equations  is  often  built 
on  understanding  the  characteristics.  In  plasticity,  the  char¬ 
acteristics  can  either  be  lines  in  the  directions  of  critical  shear 
stress  (stress  formulation)  or  slip  lines  (kinematic  formula¬ 
tion).  For  plasticity  based  on  an  associated  flow  law,  these 
characteristics  are  the  same,  although  Spencer  [1]  describes  a 
“double-shearing”  model  in  which  these  two  characteristics 
coincide  for  the  case  of  non-associated  flow.  The  characteris¬ 
tic  lines  in  [2-5]  are  interpreted  as  “stress  chains”,  the  direc¬ 
tions  of  which  are  determined  by  a  local  fabric  tensor  that 
arises  naturally  from  their  analysis.  If  the  definition  of  stress 
chain  is  extended  to  mean  a  characteristic  trajectory  along 
which  a  particular  stress  state  is  computed,  the  characteris¬ 
tics  of  limit  plasticity  and  of  media  in  the  isostatic  state  (e.g. 
[2-5])  can  both  be  thought  of  as  stress  chains.  Indeed  if  the 
boundary  conditions  causing  stress  chains  in  the  analysis  by 
[2-5]  are  applied  to  a  medium  at  the  plastic  limit  state,  stress 
chains  are  likewise  produced.  Therefore,  stress  chains  arise 
because  the  continuum  equations  are  of  hyperbolic  type. 

To  what  extent  then  are  stress  chains  similar  to  the  force 
chains  observed  in  discrete  element  simulations  of  granular 
materials  (e.g.  [7])  and  experiments  with  photoelastic  disks 


(see  [8,9])?  Here  the  relationship  between  the  directions  of 
the  characteristics  and  those  of  the  principal  stress  warrants 
attention.  Consider  the  biaxial  test  wherein  the  principal  axes 
remain  parallel  to  the  sides  of  the  specimen.  Characteristics 
that  define  stress  chains  propagate  at  oblique  angles  relative 
to  the  specimen  sides  to  straddle  a  “cone  of  influence”  [5]. 
Yet,  as  shown  in  [5],  the  principal  axes  align  locally  with 
the  direction  of  the  stress  chains.  The  situation  is  illustrated 
in  Fig.  1.  Here,  a  boundary  load  /,  acting  on  a  boundary 
particle,  results  in  the  formation  of  major  load-bearing  par¬ 
ticles  arranged  in  a  quasi-linear  fashion,  i.e.  the  force  chains 
(dark  shaded  particles  in  Fig.  1).  The  lines  correspond  to  the 
stress  chains  or  characteristics  computed  from  the  contin¬ 
uum  equations.  Similar  chains  can  be  conceived  emanating 
from  all  boundary  forces,  from  which  the  average  stress  is 
computed.  The  top  right  inset  in  Fig.  1  shows  the  stress  prop¬ 
agated  along  the  stress  chain,  where  it  is  seen  that  the  local 
principal  stress  axis  is  aligned  with  the  stress  chain  direc¬ 
tion.  This  correspondence  is  consistent  with  the  algorithm 
used  to  determine  if  a  particle  belongs  to  a  particular  force 
chain  (e.g.  [10] ).  On  the  other  hand,  the  principal  axes  of  the 
average  stress  state,  determined  by  integrating  the  boundary 
forces,  is  aligned  with  the  specimen  boundaries.  As  such, 
the  principal  axes  determined  from  the  average  stress  need 
not  correspond  to  the  local  principal  stress  axes  observed  by 
following  a  particular  stress  chain. 

In  a  similar  vein,  an  important  difference  between  stress 
chains  and  force  chains  is  that  a  stress  chain  emanates  from 
a  boundary  load  along  which  the  stress  in  the  interior  caused 
by  that  load  can  be  determined.  However,  multiple  stress 
chains  can  pass  through  an  interior  point  giving  rise  to  a  stress 
state  that  is  a  combination  of  all  intersecting  stress  chains. 
In  contrast,  the  average  stresses  within  particles  comprising 
a  force  chain  are  the  result  of  the  total  loads  on  the  system; 
strictly  speaking  the  stresses  within  the  particles  cannot  be 
propagated  from  a  particular  force,  either  on  the  boundary  or 
within  the  interior. 

To  resolve  the  issue  of  whether  stress/force  chains  follow 
the  material  fabric  or,  as  observed  in  experiments,  align  with 
the  principal  stress  axis,  it  is  useful  to  consider  what  controls 
fabric.  Like  other  fabric  tensors  used  in  granular  mechan¬ 
ics,  the  fabric  tensor  of  [2-5]  is  a  measure  of  the  relative 
position  of  grains  and  arises  naturally  from  their  analysis 
of  isostatic  states.  In  the  analysis,  the  particle  assemblage  is 
partitioned  into  particle  groups  that  surround  voids.  How¬ 
ever  there  are  differences  between  the  fabric  tensor  in  [2-5] 
and  those  typically  found  in  the  micromechanics  literature. 
In  [2-5],  the  balance  of  forces  and  moments  gives  rise  to  a 
closure  equation  expressed  in  terms  of  stresses  and  elements 
of  the  fabric  tensor.  By  contrast,  the  appellation  fabric  tensor 
is  often  applied  to  the  symmetric  tensor  created  by  averag¬ 
ing  the  outer  products  of  the  contact  normals.  This  definition 
implies  a  central  reference  particle  and  its  nearest  neighbors 
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and  is  for  the  purposes  of  describing  the  extent  to  which  this 
arrangement  is  isotropic.  The  change  in  either  measure  of 
fabric  implies  a  change  in  the  contact  topology  which,  in 
turn,  leads  to  irreversible  (plastic)  strains. 

Experiments  show  that  under  monotonic  loading  in  the 
biaxial  test,  the  principal  axes  of  the  fabric  tensor  align  with 
the  principal  axes  of  stress  as  do  the  force  chain  directions 
[7,9,11].  If  force  chains  and  stress  chains  are  controlled  by 
fabric,  how  is  the  observed  relationship  between  force  chains 
and  principal  stress  axes  [11]  explained?  The  fabric  tensor  is 
presumably  related  to  stress  history  rather  than  stress.  From 
experimental  soil  mechanics,  it  has  been  observed  that  large 
plastic  strains  occur  with  the  rotation  of  principal  axes  [12], 
suggesting  that  the  soil  fabric  rapidly  adjusts  to  the  prevail¬ 
ing  principal  stress  directions.  The  characteristic  lines  should 
therefore  be  related  to  the  stress  tensor  as  well  as  the  fabric 
tensor,  to  the  extent  that  the  principal  axes  of  each  coincide. 

A  previous  paper  [13]  explored  the  proposal  that  stresses 
in  granular  media  are  in  fact  carried  by  a  “strong  contact  net¬ 
work”  that  is  almost  isostatic  but  which  is  embedded  within 
the  overall  redundantly  constrained  (hyperstatic  or  statically 
indeterminate)  particle  system.  In  such  systems,  however, 
the  very  existence  of  the  “strong  contact  network”  is  heavily 
dependent  on  the  “weak  contact  network”;  the  latter  pro¬ 
vides  the  necessary  confining  support  to  enable  the  build-up 
of  force  along  particle  chains  that  then  result  in  the  forma¬ 
tion  of  the  former.  In  studies  of  initially  isotropic  systems 
in  2D  and  3D,  force  chains  have  been  observed  to  develop 
along  chains  of  particles  with  a  higher  degree  of  redundancy 
(hyperstatic) — specifically  through  supporting  contacts  in 
3-cycle  formations  [14,15]. 

In  this  paper,  we  advance  the  investigation  in  [13-15] 
and  the  Cosserat  constitutive  model  in  [16]  by  considering 
the  conditions  within  shear  bands  that  produce  a  statically 
determinate  system  of  continuum  equations,  even  though  the 
conditions  at  the  particle  scale  are  not  necessarily  isostatic. 
Consistent  with  the  relatively  high  porosities  found  in  shear 
bands  in  sand  [8,11,17,18],  idealized  assemblies  in  photo¬ 
elastic  disk  experiments  [9,11],  and  discrete  element  sim¬ 
ulations  [7,14]  also  exhibit  low  average  coordination  num¬ 
ber  per  particle  in  shear  bands  compared  to  the  rest  of  the 
sample,  although  this  quantity  is  observed  to  fluctuate  in  the 
band  both  spatially  and  temporally  [7,9,14].  The  rises  and 
falls  in  the  average  coordination  number — in  both  space  and 
time — are  due  to  the  continual  process  of  creation  of  new, 
and  collapse  by  buckling  of  old  force  chains,  as  illustrated 
in  Fig.  1  (bottom  right).  Corroborative  evidence  from  high 
resolution  experiments  on  sand  [17, 18]  suggest  that  regions 
of  high  (hyperstatic)  and  low  (isostatic/hypostatic)  coordi¬ 
nation  numbers  alternate  spatially  along  the  band,  with  the 
former  being  jamming  zones  where  force  chains  grow  and 
the  latter  being  unjamming  zones  where  force  chains  collapse 
by  buckling.  The  relative  dominance  of  these  two  compet¬ 


ing  events  of  jamming  and  unjamming  determines  which  of 
the  elastic  rises  and  unstable  drops  in  the  strain  evolution 
of  the  macroscopic  stress  occur  [19].  Thus,  in  the  so-called 
critical  state  regime  wherein  the  equations  of  limit  plasticity 
hold,  and  for  which  the  sample  deforms  in  the  presence  of 
fully  developed  shear  bands,  static  determinacy  in  the  contin¬ 
uum  stresses  results  from  both  spatial  and  temporal  averaging 
of  the  stress  of  their  particulate  counterparts.  Accordingly, 
this  study  explores  conditions  of  isostaticity  for  a  Cosserat 
medium  based  on  observations  of  stresses  in  the  critical  state 
in  a  manner  consistent  with  the  recent  analysis  of  confined 
force  chain  buckling  in  [16,20]. 

The  rest  of  the  paper  is  arranged  as  follows.  A  set  of  clo¬ 
sure  relations  for  an  isostatic  Cosserat  medium  is  presented  in 
Sect.  2  with  analytical  solutions  in  Sect.  3  for  the  special  case 
where  stress  characteristics  are  straight  lines.  A  numerical 
scheme  for  more  general  cases  with  an  example  is  shown  in 
Sect.  4.  To  the  extent  possible,  comparisons  are  made  with  the 
analysis  put  forth  in  Gerritsen  et  al.  [5].  Concluding  remarks 
are  given  in  Sect.  5. 

2  Closures  for  isostatic  Cosserat  continuum 

We  explore  isostatic  conditions  achieved,  on  average,  in 
space  and  time,  through  the  collective  birth  and  death  of  force 
chains  [7, 8, 1 1, 18, 19].  We  consider  a  set  of  closure  relations 
which,  when  combined  with  the  momentum  balance  relations 
for  a  Cosserat  continuum,  lead  to  a  statically  determinate  sys¬ 
tem  of  equations.  There  are  six  unknown  stress  components: 
4  for  stress  and  2  for  couple  stress.  With  three  equilibrium 
equations,  three  relations  for  the  stresses  are  needed. 

Closure  I — relation  between  normal  stresses :  A  defining 
aspect  of  the  large  strain,  fully  developed  plastic  flow  regime, 
or  “critical  state”,  for  granular  materials  undergoing  localized 
failure — is  the  almost  steady  evolution  with  strain  of  the  shear 
stress  expressed  through  the  stress  ratio  (o\\  —  022) /(&n  + 
<722 )•  The  continual  growth  of  new  force  chains  amidst  col¬ 
lapse  by  buckling  of  old  force  chains  in  the  persistent  shear 
band  is  believed  to  be  the  mechanism  underpinning  this  trend 
[7,8, 11, 17-20].  Hence  we  propose  the  closure  relation: 

022  =  Con.  (1) 

This  proportionality  relation  is  of  the  same  form  as  that 
derived  from  the  force  chain  buckling  analysis  in  [20] . 

Closure  II — relation  between  shear  stresses :  For  a  classi¬ 
cal  continuum,  the  stresses  are  symmetric  with  <712  =  <X2i- 
We  wish  to  investigate  the  influence  of  asymmetry  in  the 
stress  tensor,  when  there  is  a  gradient  in  moments,  which 
we  assume  to  be  small.  Thus  we  propose  a  simple  propor¬ 
tionality  between  shear  stresses  for  a  second  closure  relation, 
such  that  the  classical  relation  o\ 2  =  <721  is  recovered  when 
7  =  1: 
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*2  (a)  x2  (b) 


Fig.  2  Illustration  of  the  determination  of  stresses  from  characteris¬ 
tics:  (a)  through  a  point  inside  the  half  plane  from  two  boundary  forces, 
(b)  through  the  point  of  action  of  an  applied  boundary  stress 

G\2  =  To2\.  (2) 

Closure  III — relation  for  couple  stresses :  A  similar 
proportionality  relation  as  above  is  proposed: 

Ml  =  MfJL2-  (3) 


lines,  stresses  are  constant.  We  use  solutions  to  these  cases 
to  test  the  numerical  scheme  needed  to  solve  more  complex 
boundary  value  problems  (next  section).  The  boundary  con¬ 
ditions  for  the  two  cases  below  are:  (1)  two  concentrated 
forces  act  at  points  (xf,  xf)  and  (xf,  x%)  (Fig.  2a),  and  (2) 
a  concentrated  normal  force  an  =  sn  acts  at  point  (xf,  x%  ) 
(Fig.  2b). 

In  case  1,  two  characteristic  lines  from  the  points  of 
action  of  the  concentrated  forces  at  the  boundary  intersect 
at  (x\,  x2).  The  stresses  x2)  and  a2\(x\,  *2)  can  be 

obtained  in  terms  of  the  applied  concentrated  forces  at  the 
boundary,  keeping  in  mind  these  have  normal  and  shear  com¬ 
ponents.  Along  the  w\  characteristic  from  the  point  (xf ,  xf), 
we  have 

Xcrn(xi,  X2)  +  <721  (*1,  *2)  =  Xan(xf,  *2  )  +  ^21  (xf,  *2  )• 
Along  the  w2  characteristic  from  the  point  (xf,  xff), 
Xan(xi,x2)  -cr2  i(x\,x2)  =  Aan(vf,vf)  -  <r2i(*j\  x% ). 


The  special  case  of  M  —  0  follows  directly  from  the  parti¬ 
cle  kinematics  in  the  buckling  model  in  [16,20].  Note  that 
a  detailed  account  of  this  buckling  model  is  given  in  [20], 
including  a  full  justification  of  its  assumptions  and  a  com¬ 
prehensive  comparison  of  its  predictions  with  those  observed 
from:  (i)  force  chain  evolution  in  discrete  element  simula¬ 
tions  and  (ii)  global  trends  for  stress  from  experiments  on 
sand. 

Equilibrium  equations'.  With  the  above  closure  relations, 
the  equilibrium  equations  for  the  case  of  uniform  closure 


coefficients  G,T,M  become 

3<7ii  9  <721 

+  T  =0, 

dx\  dx2 

(4) 

9  <721  da  1 1 

+  G  =0, 

dx\  dx2 

(5) 

dii2  3ijl2 

Mf—  +  =  (T  —  i)(T21. 

OX 1  OX 2 

(6) 

3  Analytical  solutions  for  Cauchy  stresses  with  straight 
characteristics 

The  governing  system  of  equations  are  (4)-(6):  in  general 
G  >  0  and,  for  T  >  0,  (4)-(5)  are  strictly  hyperbolic  and 
the  influence  of  an  asymmetric  stress  tensor  can  be  probed 
for  T  >  0  (T  7^1).  Here,  we  consider  the  special  case 
M  —  0,  in  line  with  the  force  chain  buckling  model  in  [20]. 
Equation  (6)  can  be  used  to  determine  the  couple  stress  [i2 , 
once  <72i  is  established  from  (4)-(5).  In  cases  of  constant 
coefficients  G  and  T ,  and  hence  constant  GT ,  the  character¬ 
istics  comprise  two  families  of  straight  parallel  lines,  w  1  and 
w2,  with  slopes  ±X  =  d=V GT .  Along  these  characteristics 


In  case  2,  the  stresses  are  zero  except  for  the  points  along 
the  pair  of  characteristics  passing  through  (xf,  xJf ):  the  w\ 
characteristic  is  x2  =  x^  +  X(x\  —  xf)  and  the  w2  charac¬ 
teristic  is  x2  =  X2  —  X(x\  —  xf).  At  a  point  (x\,  x2)  along 
w  1,  Xa\\(x\,  x2)  +  <721  (x\ ,  x2)  =  Xs n;  while Xan(x\,  x2)  — 
<721  (*1,  *2)  =  0  along  the  w2  characteristic,  depicted  as  a 
broken  line  in  Fig.  2b.  Reminiscent  of  the  so-called  stress 
chains  in  Gerritsen  et  al.  [4,5],  the  stresses  here  are  also  con¬ 
centrated  along  the  pair  of  characteristics  through  (xf,  xf1): 
an  =  s n/2  and  <721  =  Xs n/2  for  w\,  and  an  =  ^n/2  and 
a2\  =  —Xsn/2  for  w2. 


4  Numerical  solutions 

If  the  closure  coefficients  are  general  functions  G(x \,x2), 
T(x  \,x2)  and  M(x\ ,  x2),  the  governing  system  of  equations 
are  coupled  and  a  numerical  scheme  must  be  employed. 
For  the  constant  closure  coefficients  considered  in  (4)-(6), 
the  solution  scheme  proceeds  in  three  steps:  first,  equa¬ 
tions  (4)  and  (5)  are  solved  for  an  and  <721;  second,  a22 
and  <7 12  are  determined  from  the  first  two  closure  rela¬ 
tions;  and  third  pi2  is  obtained  from  (6).  Couple  stress  /x  1 
is  obtained  from  the  third  closure  equation  (3),  but  here 
we  focus  on  the  case  M  =  0.  Equations  (4)  and  (5)  are 
solved  numerically  by  first  approximating  the  derivatives 
with  respect  to  x2  by  finite  differences;  this  yields  a  system 
of  ordinary  differential  equations  with  respect  to  x\  which 
is  then  solved  using  the  4th-order  Runge-Kutta  method. 
We  tested  this  scheme  by  solving  the  equations  of  Ger- 
ritsen’s  et  al.  [5].  Two  examples  are  considered  for  the 
following  region  and  boundary  conditions.  The  region  lies 
in  0  <  x\  <  1,  0  <  x2  <  3  and  is  represented  by  a 
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uniform  mesh  of 200  by  600  nodes  in  the  x\  and  X2  directions, 
respectively.  The  boundary  conditions  on  x\  =  Oare:  /i\  =  0 
(from  the  closure  relation),  a2i  =  0,  and  the  stress  o\\  is  pre¬ 
scribed  to  be  o\\  =  —  20N/m2  at  node  (x\,  X2)  =  (0,  1.5), 
with  o\\  =  —4.0  N/m2  at  two  adjacent  nodes  symmetrically 
located  on  either  side  of  node  (0,  1.5),  and  o\\  =  0  for  all 
other  points. 

Example  1  The  uniform  closure  coefficients  are  G  =  3,  T  = 
1  and  M  —  0.  The  resulting  characteristic  lines  have  slopes 
a/3  and  —a/3,  enabling  a  direct  comparison  with  the  ana¬ 
lytical  solution  in  Sect.  3.  In  particular,  for  the  boundary 
condition  /x 2  =  0  on  X2  =  0,  equation  (6)  leads  to  /x2  =  0 
for  the  whole  region  which  is  consistent  with  the  classical 
continuum  case.  Further,  the  results  here  can  be  compared 
with  those  in  Gerritsen  et  al.  [5]  for  that  case  when  their 
fabric  tensor  p\\  =  1,  P22  =  —3  and  pn  =  0  (see  [5] 
for  the  notation).  The  stress  cr\\  is  concentrated  along  two 
straight  bands,  of  slopes  \/3  and  —  y/3  (similar  to  the  char¬ 
acteristic  lines  for  case  2  in  Sect.  3).  These  bands  form  an 
“influence  wedge”,  reminiscent  of  the  influence  cone  by  Ger¬ 
ritsen  [5].  Like  the  influence  cone  of  Gerritsen  et  al.  [5], 
we  find  subregions  in  the  influence  wedge  in  which  o\\  is 
nonzero  or  even  positive  (i.e  tensile),  albeit  its  magnitude  is 
small.  The  analytical  solution  in  Sect.  3,  however,  dictates 
that  the  stress  inside  this  influence  wedge  is  zero.  We  ana¬ 
lyzed  the  sensitivity  of  the  numerical  solution  to  mesh  size 
and  found  that  the  magnitude  of  the  stress  in  these  subre¬ 
gions  decreases  as  the  mesh  size  used  to  approximate  the 
partial  derivative  with  respect  to  v2  decreases,  suggesting 
that  these  tensile  stress  zones  are  an  artifact  of  the  numerical 
scheme. 

Example  2  In  the  DEM  simulations  of  [7],  the  shear  stress 
fluctuates  about  a  near  constant  value;  for  the  specific  sys¬ 
tem  in  Fig.  3  of  [7],  the  stress  ratio  is  (crn  —  a22)/(crn  + 
022 )  =  0.35,  i.e.  02 2  =  2.1am  Solutions  corresponding 
to  this  relation,  with  G  =  2.1  and  M  =  0,  for  the  values  of 
T  =  0.5,  1.0,  1.5  are  shown  in  Figs.  3  and  4:  here  the  nonzero 
stress  inside  the  influence  wedge  formed  by  the  characteris¬ 
tics  is  evident,  and  most  pronounced  for  T  =  0.5. 

Finally,  Gerritsen  et  al.  [4,5]  showed  that  consideration 
of  non-constant  closure  coefficients  lead  to  “stress  leakage”, 
i.e.  non-zero  stresses  develop  inside  the  influence  cone,  the 
magnitude  of  which  depends  on  the  specific  form  of  the  clo¬ 
sure  coefficients.  For  example,  they  found  that  periodically 
varying  closure  coefficients  lead  to  a  periodic  pattern  in  the 
stresses,  and  the  two  bands  of  concentrated  stress  bear  a  peri¬ 
odic  waveform.  We  also  observed  stress  leakage  from  consid¬ 
eration  of  a  small  variation  in  our  closure  coefficients  (data 
not  shown). 


(a)  T=B,5 


0  0.1  0,2  0*3  0,4  0.5 


(b)  T=i,o 


0  04  0*2  0,3  0,4  0,5 


(C)  T=I,5 


8  0,1  0,2  0,3  8,4  8,5 

Fig.  3  Distribution  of  stress  component  o\\  for  a  constant  stress  ratio 
of  0.35  (in  accordance  with  [7])  and  a  T  =  0.5,  b  T  =  1,  c  T  =  1.5 


(a)  t=o,5 


0  0*1  0*2  0*3  0*4  0*5 


(b)  t=i,0 


0  0*1  0,2  0*3  0,4  0*5 


(c)  1=1,0 


0  8,1  0,2  0,3  0,4  B.5 

Fig.  4  Distribution  of  stress  component  <721  for  a  constant  stress  ratio 
of  0.35  (in  accordance  with  [7])  and  a  T  =  0.5,  b  T  =  1,  c  T  =  1.5 


5  Concluding  remarks 

The  static  determinacy  of  granular  media  in  the  fully  devel¬ 
oped  plastic  or  critical  state  regime  is  investigated.  The  pre¬ 
dictions  of  the  Cosserat  model  display  many  similarities  with 
those  arising  from  the  isostatic  state  investigated  by  Blu- 
menfeld  and  co-workers  [2-5].  In  both  cases,  when  the  clo¬ 
sure  coefficients  are  constant,  the  characteristics  associated 
with  a  boundary  force  define  an  influence  wedge.  Outside  the 
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wedge,  the  stress  is  zero;  inside  the  wedge,  the  stresses  are 
non-zero  albeit  these  decrease  in  magnitude  with  decreasing 
mesh  size.  We  also  observed  what  they  referred  to  as  “stress 
leakage”  for  non-constant  closure  coefficients.  These  sim¬ 
ilarities  mostly  result  from  the  governing  equations  being 
hyperbolic  for  both  the  isostatic  state  and  the  critical  state; 
however  the  physical  processes  rendering  the  systems  hyper¬ 
bolic  are  quite  different.  In  the  isostatic  case,  the  system  is 
statically  determinate  at  the  particle  scale,  which  makes  it 
possible  to  develop  continuum  equations  that  are  statically 
determinate;  at  the  critical  state,  statical  determinacy  of  the 
continuum  equations  results  from  the  spatial-temporal  aver¬ 
aging  of  forces  associated  with  the  concomitant  birth,  and 
collapse  by  buckling,  of  force  chains.  The  forces  in  the  crit¬ 
ical  state,  at  the  particulate  scale,  are  not  statically  determi¬ 
nate.  The  governing  equations  for  the  two  cases  thus  repre¬ 
sent  different  physical  processes.  Despite  the  similarity  in 
their  characteristic  solutions,  they  differ  by  their  information 
content.  The  stress  chains  produced  by  the  isostatic  condi¬ 
tion  are  similar  to  the  force  chains  seen  in  granular  media. 
In  the  case  of  the  former,  the  stress  chain  carries  informa¬ 
tion  derived  entirely  from  the  boundary  load  from  which  it 
emanates.  In  contrast,  a  force  chain  represents  the  combined 
result  of  the  boundary  force  and  the  loading  of  all  particles 
incepting  the  chain.  The  characteristics  associated  with  the 
governing  equations  for  the  critical  state  represent  the  limit¬ 
ing  stress  state,  which  is  essentially  equivalent  to  the  case  of 
limiting  plasticity.  The  stress  at  the  critical  state  is  not  neces¬ 
sarily  constant  but  fluctuates,  with  rises  and  falls  tied  to  the 
jamming  and  unjamming  occurring  inside  the  material.  The 
limiting  state  thus  represents  a  time-averaged  quantity. 

At  the  critical  state,  the  stress  states  are  episodic,  with 
rises  in  stress  related  to  force  chain  creation  (jamming)  fol¬ 
lowed  by  unstable  falls  in  stress  related  to  force  chain  buck¬ 
ling  (unjamming).  This  inference  on  mechanisms  associated 
with  jamming-unjamming  process  comes  from  direct  obser¬ 
vations  of  the  shear  stress  through  the  stress  ratio,  extracted 
from  DEM  simulations  and  experiments  on  photoelastic  disk 
assemblies  and  on  sand  [7-9,11,14,17-20].  Likewise,  the 
closure  equations  used  to  obtain  the  governing  equations 
depend,  to  some  extent,  on  direct  observation  of  how  quanti¬ 
ties  vary  in  simulations  and  experiments.  Our  ultimate  goal  is 
to  develop  a  micromechanical  model  that  predicts  the  above- 
mentioned  observations,  based  on  closure  equations  that  are 
derived  from  a  combined:  (i)  stochastic  analysis  of  the  birth- 
death  process  of  an  assembly  of  force  chains  in  distinct  stages 
of  buckling,  and  (ii)  a  comprehensive  structural  mechanics 
model  of  the  different  stages  in  the  evolution  of  a  laterally 
confined  force  chain  under  axial  compression.  This  work  is 
underway. 

Finally,  we  close  with  questions  arising  from  this  study. 
The  comparison  between  particle  scale  mechanisms  and  the 
continuum  description  produced  for  the  isostatic  and  criti¬ 


cal  states  begs  the  question  of  how  much  of  the  particulate 
behavior  survive  the  passage  to  the  continuum  representation 
via  averaging  over  the  representative  volume.  In  the  case  of 
the  critical  state,  is  the  continuum  model  even  a  proper  rep¬ 
resentative  of  the  particulate  behavior?  Indeed,  the  averag¬ 
ing  of  the  particulate  equations  entirely  masks  the  turmoil 
that  actually  occurs  at  the  particle  scale.  Such  a  difference 
might  be  significant  if,  for  example,  coupled  processes  on  the 
mesoscale  were  under  consideration  (e.g.  soil-water  interac¬ 
tion).  At  variance  to  classical  continua,  should  the  rheology 
of  complex  media  be  formulated  at  the  level  of  basic  units  (i.e. 
particle  scale) — prior  to  averaging  to  a  continuum  approxi¬ 
mation?  If  so,  what  of  the  interactions  between  mesoscopic 
entities  made  up  of  the  fundamental  units  (i.e.  force  chains)? 
In  regard  to  the  continuum  equations  describing  the  isostatic 
state,  a  similar  critical  analysis  can  be  made.  Each  stress 
chain  represents  the  stress  caused  by  a  concentrated  load 
at  the  specimen  boundary.  The  stress  state  at  any  particular 
point  is  the  result  of  adding  the  stresses  from  the  preponder¬ 
ance  of  stress  chains  intersecting  the  point.  Therefore,  the 
stress  distribution  using  the  hyperbolic  governing  equations 
of  the  isostatic  medium  is  the  net  result  of  all  boundary  forces. 
For  example,  even  though  each  stress  chain  is  aligned  in  the 
direction  of  the  principal  axis  caused  by  the  boundary  force, 
it  need  not  be  aligned  with  the  principal  axis  of  the  total 
stress  at  the  point.  In  addition,  an  increment  of  load  causes 
rearrangement  of  particles  that  is  manifest  by  a  strain  of  the 
continuum  element.  Presumably,  the  stress  and  strain  incre¬ 
ment  can  be  related  through  a  constitutive  relationship,  which 
when  combined  with  the  equations  of  equilibrium,  result  in 
elliptic  governing  equations.  The  ability  to  relate  stress  and 
strain  does  not  depend  on  the  specimen  being  in  an  isostatic 
state  (with  rigid  particles)  because  the  strain  is  caused  by 
relative  movements  among  particles  and  not  the  deformation 
of  the  particles  themselves.  It  follows  that  the  stress  can  be 
determined  either  by  hyperbolic  equations  or  elliptic  equa¬ 
tions,  which  in  principle  should  produce  the  same  results. 
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